EXO1 : Enseigner la statistique avec R

AFROMPA’R Bouaké 2025

On se propoe de montrer l’intérêt de R pour l’enseignement de la statistique aux étudiants de licence de géographie… Cet exercice correspond typiquement au travail qu’on pourrait donner à des étudiants de licence 2 ou 3 ayant suivi des cours de statistique univariée et bivariée et ayant des connaissances de base en cartographie thématique.

On propose de réaliser une analyse de la distribution de l’Indicateur de Développement humain des régions de Côte d’Ivoire entre 2002 et 2022 à l’aide des données accessibles sur le site Global Data Lab . Le programme est rédigé de telle sorte qu’on puisse facilement refaire l’exercice en prenant d’autres dates (ex. Evolution entre 2012 et 2022) ou d’autres pays d’Afrique de l’Ouest.

Code

library(sf)
library(dplyr)

(1) PREPARATION DES DONNEES STATISTIQUES

• Consigne : Après avoir chargé les fichiers idh2002.csv et idh2022.csv, construisez un tableau de donnée comportant le code et le nom des régions ainsi que leur population et leur indicateur de développement humain aux deux dates.

Code

# Importe les données
tab1<-read.table("data/exo_stats/idh1992.csv", 
                dec=".",
                sep = ";",
                header=TRUE)
tab2<-read.table("data/exo_stats/idh2022.csv", 
                dec=".", 
                sep=";",
                header=TRUE)

# Sélectionne les variables
tab1 <- tab1 %>% select(code=code, pays=pays, nom=nom, pop1992 = pop, idh1992 = idh)
tab2 <- tab2 %>% select(code=code,  pop2022 = pop, idh2022 = idh)

# Fusionne les deux tableaux
don <- left_join(tab1,tab2)


# Affiche
library(knitr)
kable(don, 
      main = "Tableau de départ")

code	pays	nom	pop1992	idh1992	pop2022	idh2022
CIVr101	CIV	Abidjan	2728	0.564	5984	0.569
CIVr102	CIV	Yamoussoukro	112	0.623	245	0.619
CIVr103	CIV	Bas Sassandra	529	0.406	1160	0.520
CIVr104	CIV	Comoe	749	0.378	1642	0.582
CIVr105	CIV	Denguele	381	0.332	836	0.440
CIVr106	CIV	Goh-Djiboua	686	0.424	1504	0.551
CIVr107	CIV	Lacs	1098	0.389	2409	0.542
CIVr108	CIV	Lagunes	1023	0.456	2245	0.566
CIVr109	CIV	Montagnes	1365	0.380	2994	0.492
CIVr110	CIV	Sassandra-Marahoue	1614	0.407	3541	0.512
CIVr111	CIV	Savanes	628	0.409	1376	0.471
CIVr112	CIV	Vallee du Bandama	520	0.443	1140	0.554
CIVr113	CIV	Woroba	574	0.199	1259	0.423
CIVr114	CIV	Zanzan	834	0.320	1829	0.491

(2) STATISTIQUE UNIVARIEE

Paramètres principaux

• Consigne : Etudiez l’évolution de l’IDH entre 1992 et 2022 en vous servant de paramètres principaux (valeurs centrales, paramètres de dispersion). Puis établissez deux histogrammes permettant de visualiser l’évolution.

Code

# sélectionne les variables
sel <- don[,c("idh1992","idh2022")]

# Tableau standard
quant<-apply(sel,2,quantile)
moy<-apply(sel,2,mean)
ect<-apply(sel,2,sd)
cv<-100*ect/moy
tab<-rbind(quant,moy,ect,cv)
row.names(tab) <-c("Minimum","Q1","Médiane","Q3","Maximum","Moyenne","Ecart-type", "C.V. (%)")

kable(tab, caption="Paramètres principaux", digits =2,
      col.names = c("Situation en 1992","Situation en 2022"),
      )

Paramètres principaux

	Situation en 1992	Situation en 2022
Minimum	0.20	0.42
Q1	0.38	0.49
Médiane	0.41	0.53
Q3	0.44	0.56
Maximum	0.62	0.62
Moyenne	0.41	0.52
Ecart-type	0.10	0.06
C.V. (%)	24.77	10.65

Histogrammes

• Consigne : Etablissez deux histogrammes permettant de visualiser la forme de la distribution de l’IDH en 1992 et en 2022

Code

par(mfrow=c(2,1))

mintot <-min(c(sel$idh1992, sel$idh2022))
maxtot <-max(c(sel$idh1992, sel$idh2022))

# Histogramme
hist(sel$idh1992,
     breaks=quantile(sel$idh1992),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 1992",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$idh1992, col="black", lwd=2)
lines(density(sel$idh1992), lty=3,lwd=2)

hist(sel$idh2022,
     breaks=quantile(sel$idh2022),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2022",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$idh2022, col="black", lwd=2)
lines(density(sel$idh2022), lty=3,lwd=2)

(3) PREPARATION DES DONNEES GEOMETRIQUES

• Consigne : Après avoir chargé le fonds de carte affichez le contour des régions et le code des unités.

Code

# Importe le fonds de carte
map <- st_read("data/exo_stats/mapreg.shp")

Reading layer `mapreg' from data source 
  `/Users/claudegrasland1/worldregio/afromapr/data/exo_stats/mapreg.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 14 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -8.599302 ymin: 4.361807 xmax: -2.494888 ymax: 10.73664
Geodetic CRS:  WGS 84

Code

# Projette
map <- st_transform(map, 2042)


# Affichage du fonds de carte
par(mar=c(0,0,3,0))
plot(map$geometry, 
     col="gray90",
     main = "Code des unités spatiales de la zone d'étude")

# Ajout du code des unités spatiales
coo<-st_coordinates(st_centroid(map))
text(coo, map$code, cex=0.5,col="black",)

(4) CARTOGRAPHIE THEMATIQUE

Cartes de stock

• Consigne : Réalisez deux cartes de stock décrivant le nombre d’habitants en 2002 et 2022 Vous utiliserez la même échelle de taille pour rendre les deux cartes comparables.

Code

library(mapsf)
map<-map[,c("code","geometry")]
map_don <- merge(map, don, by="code")

maxequ<-max(don$pop012,don$pop2022)

par(mfrow=c(1,2))
mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop1992",
       val_max = maxequ, inches=0.1, col="gray20", 
       leg_title = "Population (milliers)",)
mf_layout(title="1992",frame = T, credits = "Source : Global Data Lab, 2025")

mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop2022",
       val_max = maxequ, inches=0.1, col="gray20",
       leg_title = "Population (milliers)")
mf_layout(title="2022",frame = T, credits = "Source : Gloabl Data Lab 2025")

Cartes de ratio (choroplèthes)

• Consigne : Réalisez deux cartes de taux décrivant le niveau de l’IDH en 2002 et 2022. Pour les rendre comparables vous utiliserez dans chaque carte une partition en quartiles (4 classes d’effectifs égaux)

Code

library(mapsf)
map_don <- merge(map, don, by="code")
maxequ<-max(don$hdi2002,don$hdi2022)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="idh1992",
       breaks = "quantile",nbreaks = 4, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="1992",frame = T, credits = "Source : Global Data Lab, 2025")

mf_map(map_don, type="choro", var="idh2022",
       breaks = "quantile",nbreaks = 4, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="2022",frame = T, credits = "Source : Global Data Lab, 2025")

(5) STATISTIQUES BIVARIEES

Nuage de points

• Consigne : Tracez un nuage de point montrant l’évolution de l’indicateur de développement humain entre les deux dates.

Code

# prépration de l'analyse
code<-don$code
nom<-don$nom
X<-don$idh1992
Y<-don$idh2022
tab<-data.frame(code,nom,X,Y)

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 1992",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

Analyse de la corrélation

• Consigne : calculez les coefficients de corrélation de Pearson et Spearman et testez leur sgnificativité.

Code

cor.test(X,Y, method="pearson")

    Pearson's product-moment correlation

data:  X and Y
t = 4.9223, df = 12, p-value = 0.0003524
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5073463 0.9403486
sample estimates:
      cor 
0.8177876 

Code

cor.test(X,Y, method="spearman")

    Spearman's rank correlation rho

data:  X and Y
S = 134, p-value = 0.006405
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7054945 

Droite de régression

• Consigne : calculez l’equation de la droite de régression et tracez- là sur le graphique.

Code

modreg <- lm(Y~X)
summary(modreg)

Call:
lm(formula = Y ~ X)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.052586 -0.016570 -0.001574  0.019274  0.072368 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.33948    0.03848   8.823 1.36e-06 ***
X            0.45013    0.09145   4.922 0.000352 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03342 on 12 degrees of freedom
Multiple R-squared:  0.6688,    Adjusted R-squared:  0.6412 
F-statistic: 24.23 on 1 and 12 DF,  p-value: 0.0003524

Code

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 1992",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

abline(modreg,col="black",lwd=1)

Analyse des résidus

• Consigne : Calculez les valeurs théoriques prévus par le modèle de régression et les résidus. Affichez le tableau correspondant après l’avoir trié par ordre de résidus croissants.

Code

tab$Y_est <- modreg$fitted.values
tab$Y_res <- modreg$residuals
tab<-tab[order(tab$Y_res),]
kable(tab, digits=3)

	code	nom	X	Y	Y_est	Y_res
11	CIVr111	Savanes	0.409	0.471	0.524	-0.053
5	CIVr105	Denguele	0.332	0.440	0.489	-0.049
1	CIVr101	Abidjan	0.564	0.569	0.593	-0.024
9	CIVr109	Montagnes	0.380	0.492	0.511	-0.019
10	CIVr110	Sassandra-Marahoue	0.407	0.512	0.523	-0.011
13	CIVr113	Woroba	0.199	0.423	0.429	-0.006
3	CIVr103	Bas Sassandra	0.406	0.520	0.522	-0.002
2	CIVr102	Yamoussoukro	0.623	0.619	0.620	-0.001
14	CIVr114	Zanzan	0.320	0.491	0.484	0.007
12	CIVr112	Vallee du Bandama	0.443	0.554	0.539	0.015
6	CIVr106	Goh-Djiboua	0.424	0.551	0.530	0.021
8	CIVr108	Lagunes	0.456	0.566	0.545	0.021
7	CIVr107	Lacs	0.389	0.542	0.515	0.027
4	CIVr104	Comoe	0.378	0.582	0.510	0.072

Cartographie des résidus

• Consigne : Cartographiez les résidus après les avoir standardisés.

Code

library(mapsf)

# Standardisation des résidus
tab$Y_res_std<-tab$Y_res/sd(tab$Y_res)

# Jointure avec la carte
map<-map[,c("code","geometry")]
map_reg <- merge(map, tab, by="code")

# Choix de la palette et des classes
library(RColorBrewer)
mypal<-brewer.pal(n = 6, name = "RdYlBu")
mybreaks = c(-10, -2,-1,0,1,2,10)

mf_map(map_reg, type="choro", var="Y_res_std",
       pal = mypal, breaks=mybreaks,
       leg_title = "Résidus standardisés",leg_val_rnd = 1)
mf_layout(title="Ecarts à la tendance 1992-2022",frame = T, credits = "Source : Global Data Lab, 2025")

(6) ANALYSE SPATIALE

On se propose d’analyser l’effet de la distance à la capitale économique (Abidjan) ou politique (Yamoussoukro) sur l’évolution de lIDH entre 1992 et 2022

Distance aux deux capitales.

• Consigne : ajoutez au tableau de données une colonne correspondant à la distance en km à Abidjan ou Yamoussoukro et faites en une cartographie en prenant comme bornes de classes 0, 100, 200, 300, 400, 500, 1000 km.

Code

map_don$capeco<-0
map_don$capeco[1]<-1
capeco<-map_don[map_don$capeco==1,]
map_don$distcapeco<-1+as.numeric(st_distance(st_centroid(map_don),st_centroid(capeco)))/1000

map_don$cappol<-0
map_don$cappol[2]<-1
cappol<-map_don[map_don$cappol==1,]
map_don$distcappol<-1+as.numeric(st_distance(st_centroid(map_don),st_centroid(cappol)))/1000



# Choix de la palette et des classes
mypal<-brewer.pal(n = 7, name = "Greys")
mybreaks = c(0,100,200,300,400,500, 1000)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="distcapeco",
       pal = mypal, breaks=mybreaks,
       leg_title = "en km",leg_val_rnd = 0)
mf_layout(title="Distance à la capitale économique",frame = T, credits = "Source : EE CIST 2025")

mf_map(map_don, type="choro", var="distcappol",
       pal = mypal, breaks=mybreaks,
       leg_title = "en km",leg_val_rnd = 0)
mf_layout(title="Distance à la capitale politique",frame = T, credits = "Source : EE CIST 2025")

Relation entre développement et distance à la capitale

• Consigne : Etudiez les corrélations entre l’IDH en 1992 (Y1) ou 2022 (Y2) et la distance à la capitale économique (X1) ou politique (X2). Que pouvez-vous en conclure ?

Code

X1<-map_don$distcapeco
X2<-map_don$distcappol
Y1<-map_don$idh1992
Y2<-map_don$idh2022

par(mfrow=c(2,2))
plot(X1,Y1,main="Distance à Abidjan 1992", sub = round(cor(X1,Y1),3), pch=20)
abline(lm(Y1~X1),col="red")
plot(X2,Y1, main = "Distance à Yamoussoukro 1992",sub = round(cor(X2,Y1),3),pch=20)
abline(lm(Y1~X2),col="red")
plot(X1,Y2,main="Distance à Abidjan 2022", sub = round(cor(X1,Y2),3), pch=20)
abline(lm(Y2~X1),col="red")
plot(X2,Y2, main = "Distance à Yamoussoukro 2022",sub = round(cor(X2,Y2),3),pch=20)
abline(lm(Y2~X2),col="red")